#required packages to make a packagepacman::p_load(devtools, roxygen2, testthat)#dependecies#pacman::p_load(tidyverse, data.table, DT, plotly, gridExtra, latex2exp, here, logr, future.apply, furrr)#pacman::p_load(dplyr,parallel,sandwich,here,future,logr,profvis,lobstr,parallel)#https://ourcodingclub.github.io/tutorials/writing-r-package/#create_package#https://combine-australia.github.io/r-pkg-dev/functions.html#load_all()#run these two lines and then library(GridMaker) in normal rstudio to install package# load_all loads a package. It roughly simulates what happens when a package is installed and loaded with library().devtools::load_all()devtools::install()
Code
library(vars)
Loading required package: MASS
Loading required package: strucchange
Loading required package: zoo
Attaching package: 'zoo'
The following objects are masked from 'package:base':
as.Date, as.Date.numeric
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ stringr::boundary() masks strucchange::boundary()
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
✖ dplyr::select() masks MASS::select()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
Code
library(plotly)
Attaching package: 'plotly'
The following object is masked from 'package:ggplot2':
last_plot
The following object is masked from 'package:MASS':
select
The following object is masked from 'package:stats':
filter
The following object is masked from 'package:graphics':
layout
Code
library(iStabi)
Warning: replacing previous import 'purrr::transpose' by
'data.table::transpose' when loading 'iStabi'
Warning: replacing previous import 'data.table::first' by 'dplyr::first' when
loading 'iStabi'
Warning: replacing previous import 'data.table::last' by 'dplyr::last' when
loading 'iStabi'
Warning: replacing previous import 'data.table::between' by 'dplyr::between'
when loading 'iStabi'
varModel <-VAR(cbind('i'= sam$y1, 's'= sam$y2), p =1, type ="const")Res =residuals(varModel)
Code
# Use for computing the moment functionOmega <-cbind('ii'= Res[, 'i']*Res[, 'i'], 'is'= Res[, 'i']*Res[, 's'],'ss'= Res[, 's']*Res[, 's'])Y <-cbind( 'omega_ii'= Omega[, 'ii'], 'omega_ss'= Omega[,'ss'], 'omega_is'= Omega[,'is'] )# Compute the moment function at each t. # INPUT: a named vector with values for the structural parameters# OUTPUT: a TxG matrix where each row corresponds to the value of the moment function on that obseration (date)# an each column the value fot ehfat moment equationGET_MOM_VALUES <-function(t0, ...) {list2env(as.list(t0),envir =environment())#Y <- cbind( 'omega_ii' = Res[,'i']^2, 'omega_ss' = Res[,'s']^2, 'omega_is' = Res[,'i'] * Res[,'s'] )#Y %*% c(-theta0['alpha'],# -theta0['beta'],# 1 + theta0['alpha'] * theta0['beta'])#fT as.matrix(-alpha*Omega[, 'ii'] - beta*Omega[, 'ss'] + (1+alpha*beta) * (Omega[, 'is']))# fT - lambda#matrix(fT - mean(fT),# ncol = 1) # (T-varLeng) X 1}
#### Search settings ----NCORES = params$NCORESSTEPSPERCORE =200# Size of adjacent cells to explores is NCORES * STEPSPERCORE# PARAMS_CFG a tibble where each row defines properties of the parameters in the name column# stepSize: size of the step for each parameter. You can have different sizes (and precision) for each paramter.# minVal and maxVal: Limits of the parameter space for each parameters. A list where each element defines the minimum and maximum for each parameter# ORDER: the order of the parameter has to match the order in theta (as used by the get_b() fcn)PARAMS_CFG <-tribble(~name, ~iniVal, ~stepSize, ~minVal, ~maxVal,'alpha', 0.3, 0.05, -1, 1,'beta', -0.2, 0.05, -1, 1) %>%mutate('N_TICKS'= ((maxVal-minVal)/stepSize)+1)
Code
NEQS =1MODELEQS <-list(list('Y'= Y, 'get_b'=NULL))#'X'=X, 'Z'=Z))NCORES =4# Level of condfidence of the S set. The closer to 1, the more robust the search is to local minima but it takes longer.SPVALUE =0.95# R is NEQS*ncol(Z) where ncol(Z) is the number of instruments# By default this matrix has to be automatically built depending on the dimensions of MODELEQS.# The user can alternatively provide this matrixR =1if (exists('Z', mode='numeric')) { R <-rbind(cbind(matrix(1, nrow=ncol(Z), ncol=ncol(Z))) )}theta0 = PARAMS_CFG$iniValnames(theta0) = PARAMS_CFG$namehead(GET_MOM_VALUES(theta0))
#logFileName <- 'log_file.txt'#log_open(logFileName, logdir=FALSE, compact=TRUE)#log_print(paste0("Local search? ", params$LOCALSEARCH))#log_print(paste0("IniRows size (in bytes): ", lobstr::obj_size(IniRows)))# rover_Perseverance will update the Grid in the Global env#plan(multisession, workers=params$NCORES)#log_print("Calling garbage collector before Stage2_rover_Perseverance")#expResult <- Stage2_rover_Perseverance(iniVal=PARAMS_CFG$iniVal, WTLags = 1, localSearch = FALSE) # updates the Grid by in placeGrid <-readRDS( 'Grid_Unidentified_model.rds')# saveRDS(Grid, 'Grid_Unidentified_model.rds')NFIXPARAMS=1
Code
head(Grid)
Key: <alpha, beta>
alpha beta S qLLStab ROW_NUMBER c2_hat genS_qLL
<num> <num> <num> <num> <int> <num> <num>
1: -1 -1.00 254.7320 11.16862 1 NA 242.7432
2: -1 -0.95 249.3581 11.29589 2 NA 237.9850
3: -1 -0.90 242.5221 11.36426 3 NA 231.8389
4: -1 -0.85 234.3027 11.37068 4 NA 224.3731
5: -1 -0.80 224.8306 11.31425 5 NA 215.7057
6: -1 -0.75 214.2824 11.19630 6 NA 205.9985
Code
Grid[order(qLLStab)][1:10]
alpha beta S qLLStab ROW_NUMBER c2_hat genS_qLL
<num> <num> <num> <num> <int> <num> <num>
1: 0.05 -1.00 22.836937 0.9207345 841 NA 21.681586
2: 0.00 -1.00 29.695267 0.9310352 801 NA 27.926732
3: 0.10 -1.00 17.258594 0.9463928 881 NA 16.636024
4: -0.05 -1.00 38.028594 0.9873345 761 NA 35.558784
5: 0.15 -1.00 12.775013 0.9995450 921 NA 12.613193
6: -0.05 -0.95 22.475247 1.0509994 762 NA 21.483042
7: -0.10 -0.95 30.765062 1.0541015 722 NA 29.022340
8: 0.20 -1.00 9.217574 1.0732271 961 NA 9.452839
9: 0.00 -0.95 15.953995 1.0996115 802 NA 15.603243
10: -0.10 -1.00 48.029135 1.1011542 721 NA 44.764004
BUILD SETS
Critical values
Critical values for \(qLL-\tilde{S}\) for 1 moment condition from table VII Suplementary material, p3 Magnusson Mavroeidis (2014).
THe S-test follows a Chi-Square distribution with 1 degrees of freedom (equatl to the number of moment conditions).
Grid has 1600 rows.
The number of fixed identified parameters is 1.
The number of estimated parameters under the null is 0.
Code
NSIPARAMS =0# not needed in the new version of the packagePCONFS =c(0.85, 0.90)get_CVs <-function() { CVs <-list('qLLStab'=NA, 'S'=NA, 'genS_qLL'=NA) CVs$qLLStab <-filter(iStabi::data_CVs$qLLStab, DG_F == KZ & P_CONF %in% PCONFS) CVs$genS_qLL <-filter(iStabi::data_CVs$genS_qLL, MOMENT_CONDITIONS == KZ & STRONGLY_IDENTIFIED_PARAMETERS == NSIPARAMS & P_CONF %in% PCONFS) CVs$S <-tibble('P_CONF'= PCONFS,'CUTOFF_TOP'=map_dbl(PCONFS, ~qchisq(.x, KZMKX))) CVs}#CV$S <- tibble('SIG' = c(0.1, 0.05, 0.01, 0),# 'VALUE' = c(0, qchisq(.9, KZMKX), # qchisq(.95, KZMKX), # qchisq(.99, KZMKX)))CVs <-get_CVs()
Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
Warning in geom_point(aes(x = -0.75, y = 0.15), size = 3, colour = "black", : All aesthetics have length 1, but the data has 546 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing
a single row.
Code
pGplot
Warning in geom_point(aes(x = -0.75, y = 0.15), size = 3, colour = "black", : All aesthetics have length 1, but the data has 546 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing
a single row.
Source Code
---title: "Unidentified case"author: "Emi Carlevaro"date: todaydate-meta: 14/01/2026editor: mode: sourceformat: html: code-fold: true code-tools: true self-contained: true page-layout: full toc: trueparams: IDSAMPLE: "" NCORES: 6 LOCALSEARCH: 'true'---## INIT```{r, eval=FALSE}#required packages to make a packagepacman::p_load(devtools, roxygen2, testthat)#dependecies#pacman::p_load(tidyverse, data.table, DT, plotly, gridExtra, latex2exp, here, logr, future.apply, furrr)#pacman::p_load(dplyr,parallel,sandwich,here,future,logr,profvis,lobstr,parallel)#https://ourcodingclub.github.io/tutorials/writing-r-package/#create_package#https://combine-australia.github.io/r-pkg-dev/functions.html#load_all()#run these two lines and then library(GridMaker) in normal rstudio to install package# load_all loads a package. It roughly simulates what happens when a package is installed and loaded with library().devtools::load_all()devtools::install()``````{r}library(vars)library(tidyverse)library(plotly)library(iStabi)``````{r}sam <- iStabi::data_sim_biVAR```## THE MODELTrue parameters are$$ \begin{align}y_{1,t} &= \textcolor{blue}{\beta} \,\, y_{2,t} + 0.5 y_{1, t-1} + 0.2 y_{2, t-1} + w_{1, t} \\y_{2,t} &= \textcolor{blue}{\alpha} \,\, y_{1,t} + 0.1 y_{1, t-1} + 0.4 y_{2, t-1} + w_{2, t} \end{align}$$$$ \begin{bmatrix}1 & -0.15 \\ 0.75 & 1\end{bmatrix} y_t = \begin{bmatrix}0.5 & 0.2 \\ 0.1 & 0.4\end{bmatrix} y_{t-1} + w_t$$with $\beta = 0.15$ and $\alpha = -0.75$.## ESTIMATION VAR```{r}varModel <-VAR(cbind('i'= sam$y1, 's'= sam$y2), p =1, type ="const")Res =residuals(varModel)``````{r}# Use for computing the moment functionOmega <-cbind('ii'= Res[, 'i']*Res[, 'i'], 'is'= Res[, 'i']*Res[, 's'],'ss'= Res[, 's']*Res[, 's'])Y <-cbind( 'omega_ii'= Omega[, 'ii'], 'omega_ss'= Omega[,'ss'], 'omega_is'= Omega[,'is'] )# Compute the moment function at each t. # INPUT: a named vector with values for the structural parameters# OUTPUT: a TxG matrix where each row corresponds to the value of the moment function on that obseration (date)# an each column the value fot ehfat moment equationGET_MOM_VALUES <-function(t0, ...) {list2env(as.list(t0),envir =environment())#Y <- cbind( 'omega_ii' = Res[,'i']^2, 'omega_ss' = Res[,'s']^2, 'omega_is' = Res[,'i'] * Res[,'s'] )#Y %*% c(-theta0['alpha'],# -theta0['beta'],# 1 + theta0['alpha'] * theta0['beta'])#fT as.matrix(-alpha*Omega[, 'ii'] - beta*Omega[, 'ss'] + (1+alpha*beta) * (Omega[, 'is']))# fT - lambda#matrix(fT - mean(fT),# ncol = 1) # (T-varLeng) X 1} ``````{r, eval=FALSE}GET_B <- function(theta0) { c(-theta0['alpha'], -theta0['beta'], 1 + theta0['alpha'] * theta0['beta'])}```## SEARCH### Search settings```{r}#### Search settings ----NCORES = params$NCORESSTEPSPERCORE =200# Size of adjacent cells to explores is NCORES * STEPSPERCORE# PARAMS_CFG a tibble where each row defines properties of the parameters in the name column# stepSize: size of the step for each parameter. You can have different sizes (and precision) for each paramter.# minVal and maxVal: Limits of the parameter space for each parameters. A list where each element defines the minimum and maximum for each parameter# ORDER: the order of the parameter has to match the order in theta (as used by the get_b() fcn)PARAMS_CFG <-tribble(~name, ~iniVal, ~stepSize, ~minVal, ~maxVal,'alpha', 0.3, 0.05, -1, 1,'beta', -0.2, 0.05, -1, 1) %>%mutate('N_TICKS'= ((maxVal-minVal)/stepSize)+1)``````{r}NEQS =1MODELEQS <-list(list('Y'= Y, 'get_b'=NULL))#'X'=X, 'Z'=Z))NCORES =4# Level of condfidence of the S set. The closer to 1, the more robust the search is to local minima but it takes longer.SPVALUE =0.95# R is NEQS*ncol(Z) where ncol(Z) is the number of instruments# By default this matrix has to be automatically built depending on the dimensions of MODELEQS.# The user can alternatively provide this matrixR =1if (exists('Z', mode='numeric')) { R <-rbind(cbind(matrix(1, nrow=ncol(Z), ncol=ncol(Z))) )}theta0 = PARAMS_CFG$iniValnames(theta0) = PARAMS_CFG$namehead(GET_MOM_VALUES(theta0))``````{r}verify_global(PARAMS_CFG, Y)Times =NROW(Res) # time periodsgenerate_global(PARAMS_CFG, SPVALUE)```### Global search grid```{r}gridLines <-build_lines(PARAMS_CFG, 'iniVal')# paramsCfg <- PARAMS_CFG# thisGridLines <- thisGridLines# NEQS = 1# psUnderNull = 1Grid <-build_grid(PARAMS_CFG, gridLines)```Test```{r}theta0 =c('alpha'=-0.75, 'beta'=0.15) # the true theta with the estimated (nuisance) parameter 0.378comp_tests(theta0, typeVarF='NeweyWest', VffLags =1)comp_tests(theta0, typeVarF='NeweyWest', VffLags =NULL)comp_tests(theta0, typeVarF='const', VffLags =NULL)comp_tests(theta0, typeVarF='HAC', VffLags =NULL)```### Start search```{r}#logFileName <- 'log_file.txt'#log_open(logFileName, logdir=FALSE, compact=TRUE)#log_print(paste0("Local search? ", params$LOCALSEARCH))#log_print(paste0("IniRows size (in bytes): ", lobstr::obj_size(IniRows)))# rover_Perseverance will update the Grid in the Global env#plan(multisession, workers=params$NCORES)#log_print("Calling garbage collector before Stage2_rover_Perseverance")#expResult <- Stage2_rover_Perseverance(iniVal=PARAMS_CFG$iniVal, WTLags = 1, localSearch = FALSE) # updates the Grid by in placeGrid <-readRDS( 'Grid_Unidentified_model.rds')# saveRDS(Grid, 'Grid_Unidentified_model.rds')NFIXPARAMS=1``````{r}head(Grid)``````{r}Grid[order(qLLStab)][1:10]```## BUILD SETS### Critical valuesCritical values for $qLL-\tilde{S}$ for `r KZ` moment condition from table VII Suplementary material, p3 Magnusson Mavroeidis (2014).THe S-test follows a Chi-Square distribution with **`r KZMKX`** degrees of freedom (equatl to the number of moment conditions).`Grid` has `r NROW(Grid)` rows.The number of *fixed* identified parameters is **`r NFIXPARAMS`**.The number of *estimated* parameters under the null is **`r NSIPARAMS`**.```{r CVs}NSIPARAMS = 0 # not needed in the new version of the packagePCONFS = c(0.85, 0.90)get_CVs <- function() { CVs <- list('qLLStab'=NA, 'S'=NA, 'genS_qLL'=NA) CVs$qLLStab <- filter(iStabi::data_CVs$qLLStab, DG_F == KZ & P_CONF %in% PCONFS) CVs$genS_qLL <- filter(iStabi::data_CVs$genS_qLL, MOMENT_CONDITIONS == KZ & STRONGLY_IDENTIFIED_PARAMETERS == NSIPARAMS & P_CONF %in% PCONFS) CVs$S <- tibble('P_CONF' = PCONFS, 'CUTOFF_TOP' = map_dbl(PCONFS, ~qchisq(.x, KZMKX))) CVs}#CV$S <- tibble('SIG' = c(0.1, 0.05, 0.01, 0),# 'VALUE' = c(0, qchisq(.9, KZMKX), # qchisq(.95, KZMKX), # qchisq(.99, KZMKX)))CVs <- get_CVs()```Critical values are:::: {.panel-tabset}### S```{r}CVs$S```### qLLStab```{r}CVs$qLLStab```### genS_qLL```{r}CVs$genS_qLL```:::### Data for plotting```{r}build_set <-function(setName) {# setName = 'qLLStab' thisCV <- CVs[[setName]] statMax =filter(thisCV, P_CONF == PCONFS[2])$CUTOFF_TOP thisPconf = thisCV$P_CONF thisCV = thisCV$CUTOFF_TOP# 546 rows theseCols <-c(PARAMS_CFG$name, setName) aSet = Grid[get(setName) <= statMax, ..theseCols][, sigAt := thisPconf[ 1+findInterval(get(setName), thisCV,left.open=TRUE)]] aSet}``````{r}Sets <-list('qLLStab'=build_set('qLLStab'),'S'=build_set('S'),'genS_qLL'=build_set('genS_qLL'))```## PLOTTING### Grid for plottingForce the starting values to be a multiple of the step size. This avoids problems and is clear what the true parameter space is.### PlotsTrue parameter vector $ \color{blue}{\theta} =(\alpha = -0.75, \beta = 0.15)$#### Global searchScatterplot for all searched values (useful for the global search)##### S plot```{r}# INTERACTIVE PLOTs ----#### S-plot ----pSbase <-plot_ly(data = Sets$S[sigAt<=0.90,], x=~alpha, y=~beta,color='grey') %>%layout(title ='S-set (90% confidence)',xaxis =list(title=plotly::TeX('\\alpha'), tickvals=seq(from=-1, to=1, by=0.20),range =c(-1, 1),fixedrange =TRUE),yaxis =list(title=plotly::TeX('\\beta'), tickvals=seq(from=-1, to=1, by=0.20),range =c(-1, 1),fixedrange =TRUE)) %>%add_trace(x =-0.75,y =0.15,type ="scatter", mode ="markers",marker =list(symbol ="square", size =12, color ="red" ),inherit =FALSE,showlegend =FALSE) %>%config(mathjax ="cdn")pSbase %>%add_trace(type='scatter', mode='markers', marker =list(color='blue'))```##### qLL-Stab plot```{r}# INTERACTIVE PLOTs ----#### S-plot ----pSbase <-plot_ly(data = Sets$qLLStab[sigAt<=0.90,], x=~alpha, y=~beta,color='grey')%>%layout(title ='qLLStab-set (90% confidence)',xaxis =list(title=plotly::TeX('\\alpha'), tickvals=seq(from=-1, to=1, by=0.20),range =c(-1, 1),fixedrange =TRUE),yaxis =list(title=plotly::TeX('\\beta'), tickvals=seq(from=-1, to=1, by=0.20),range =c(-1, 1),fixedrange =TRUE)) %>%add_trace(x =-0.75,y =0.15,type ="scatter", mode ="markers",marker =list(symbol ="square", size =12, color ="red" ),inherit =FALSE,showlegend =FALSE) %>%config(mathjax ="cdn")pSbase %>%add_trace(type='scatter', mode='markers', marker =list(color='blue'))``````{r}pBase <-ggplot(data=Sets$qLLStab[sigAt<=0.90,], mapping=aes(x=alpha, y=beta)) +scale_x_continuous(name=TeX(r"(\\alpha)"), limits=c(-1, 1), breaks=seq(from=-1, to=1, by=0.1)) +scale_y_continuous(name=TeX(r'($\\beta$)'), limits=c(-1, 1), breaks=seq(from=-1, to=1, by=0.1)) +theme_bw() +theme(panel.grid=element_line(colour='#999999', linetype='14'), panel.grid.minor=element_line(colour='white'), text=element_text(size=10), axis.title.x =element_text(colour='red'),axis.title.y =element_text(colour='blue', size=rel(0.9), hjust=0.9),axis.text.x =element_text(size=rel(1.25), angle=90, vjust=0.5), axis.text.y =element_text(size=rel(1.25))) +geom_vline(xintercept=0, colour ='black', linetype ='solid') +geom_hline(yintercept=0, colour ='black', linetype ='solid') +# Highlightinh the TRUE valuegeom_point(aes(x=-0.75, y=0.15), size=3, colour='black', fill='yellow')# ScatterplotpGplot <- pBase +geom_point(aes(colour=qLLStab))ggsave('Unidentified_model_qLLStab_set.png', plot=pGplot, width=6, height=5)pGplot```